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Summary The result after N steps of an implicit Runge-Kutta 
time discretization of an inhomogeneous linear parabolic differential 
equation is computed, up to accuracy e, by solving only 

/ 1 
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linear systems of equations. We derive, analyse, and numerically il- 
, lustrate this fast algorithm. 

in 
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In the method of lines, semi-discretization in space turns a linear 
parabolic differential equation into a large, stiff system of ordinary 
5— i ' differential equations 



1 Introduction 



u'(t)+Au(t) = g(t), u(0) = u , (1) 

possibly with a mass matrix multiplying the time derivative. This 
system is subsequently discretized in time, e.g., by the implicit Euler 
method with step size h, 

(I + hA)u n = u n -\ + hg(t n ), n = l,...,N. 
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The approximation un for a prescribed step number N is thus ob- 
tained by solving a sequence of N linear systems with a matrix of 
the form A + A, where A = 1/h in the implicit Euler method. For 
N steps with a higher-order, m-stage Runge-Kutta method, there 
are mN such linear systems, possibly with complex A as in the excel- 
lent Radau IIA methods. Even if fast techniques such as multi-grid 
methods are used, solving the linear systems of equations typically 
constitutes the main computational cost, in particular for problems 
in complicated spatial geometries. 

In this paper we propose an algorithm to compute the implicit 
Runge-Kutta approximation ujy at a fixed time T = Nh, up to an 
arbitrary accuracy e, by doing N Runge-Kutta steps for differential 
equations of the form y'(t) = Xy(t) + g(t), each step in parallel for 
0(log(l/e)) complex parameters A, and by solving only 

0(logiV log^) linear systems 

with matrices of the form A + A, all of which can be solved in par- 
allel. The constant in this work estimate is moderate: for a relative 
accuracy of 10 -5 and N < 10 5 time steps we need to solve less than 
100 linear systems! For large step numbers N, the number of linear 
systems is thus dramatically reduced, both in a sequential and in a 
parallel computational setting. 

The algorithm is highly efficient for computing Runge-Kutta ap- 
proximations to the solution of Q) at a relatively small number of 
selected time points or of short subintervals, but it is not useful for 
computing all values ui, . . . , wjv- 

Basic ingredients of the algorithm are the following: 

— the discrete variation-of-constants formula for the Runge-Kutta 
method; 

— the Cauchy integral representation of the approximations to the 
operator exponential; 

— the discretization of the contour integrals, using 0(log N) contours 
with 0(log(l/e)) quadrature points each; 

— the discrete semigroup property, which permits us to reinterpret 
the split sums as Runge-Kutta approximations to solutions of 
equations of the form y'(t) = Xy(t) + g(t). 

The algorithm given here is closely related to the fast convolution 
algorithms developed in |lllll3j . The error analysis for the discretized 
contour integrals follows the analysis of inverse Laplace transform 
approximations in 
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Discretized contour integrals have been used previously in several 
instances in the numerical solution of parabolic equations: for homo- 
geneous problems (g = 0) in ^IJ similarly to Talbot's method 19 for 
the inversion of the Laplace transform (s + A) -1 ^, and more recently 
for inhomogeneous problems |15|l5] using the Laplace transform of the 
inhomogeneity g or assuming special properties, in particular analyt- 
icity, of g. In contrast, the present algorithm works directly with the 
discrete values g(t) that are used in the Runge-Kutta discretization 
of No smoothness conditions for g are needed. This is because 
the algorithm approximates the discrete result of the Runge-Kutta 
method, with an error that does not depend on the smoothness of 
either the inhomogeneity or the solution. Of course, to make sense, 
the Runge-Kutta discretization of (pQ) with the considered step size 
h should be sufficiently accurate, which in turn does depend on the 
smoothness of g (see jlUj for Runge-Kutta error bounds for parabolic 
equations in terms of the data). 

About the differential equation (^Q) we assume that A is sectorial: 
there exist real constants M and a and an angle ip < 5 such that the 
resolvent is bounded by 



IKA + A)- 1 !! < -, for |arg(A-a)| < tt - p. (2) 

A — o~\ 



Here || • || is the operator norm corresponding to a vector norm, also 
denoted by || • ||. Clearly, for a symmetric positive semi-definite ma- 
trix A the bound ((2j) holds in the Euclidean norm with a = and 
M = 1/ sin if for any positive angle cp. More generally, condition 
P|) includes also non-symmetric operators such as those arising in 
convection-diffusion equations. In many situations resolvent bounds 
(J2J) in LP norms are known to be inherited from the continuous prob- 
lem by finite differences or finite elements, uniformly in the spatial 
discretization parameter (see, e.g., ^121 )• 

In Section |21 we review the discrete variation-of-constants formula 
for implicit Runge-Kutta methods, and in Section |3] we describe the 
discretization of the contour integrals for the rational approximations 
to the matrix exponential. The fast algorithm is given in Section 01 
including an extension to systems with a mass matrix. A numeri- 
cal example illustrates the performance of the algorithm in Section 
Finally, Section H3 analyses the error of the contour integral discretiza- 
tion, which is the only error source in the algorithm. 



4 



M. Lopez-Fernandez, C. Lubich, C. Palencia, and A. Schadle 



2 The discrete variation-of-constants formula 

In this preparatory section we recall the discrete variation-of-constants 
formula for implicit Runge-Kutta methods; cf., e.g., 

An implicit m-stage Runge-Kutta method applied to (JJ) yields, 
at t n = nh, an approximation u n to u(t n ), given recursively by 

m 

Vni = u n + hy~] djj (-Av nj + g(t n + Cjhj) , l<i<m, (3) 

3=1 
m 

u n +i = u n + h^bj (-Av nj + g(t n + Cj/t)) . (4) 

3=1 

The method is determined by its coefficients a»j, bj,Ci (i, j = 1, . . . , m). 
We denote the Runge-Kutta matrix by (2 = (aij) and the row vector 
of the weights by b T = (bj ) . Eliminating the internal stages v n i results 
in 

m 

u n+ i=r(-hA)u n + h S ^q i (-hA)g(t n + c i h), n > 0, (5) 
i=l 

where the rational approximation r(z) to e z is defined by 

r(z) = 1 + zb T (I - zQyH (6) 

with 11 = (1, . . . , 1) T , and where the rational functions qi(z) are the 
entries of the row vector 1 

q(z) = ( Ql (z), q m (z)) = b T {I - zO)-\ (7) 

We assume that the eigenvalues of the Runge-Kutta matrix (2 have 
positive real part, and that the method is L-stable, i.e., 

\r(z)\ < 1 for Rez<0, and r(oo) = 0. (8) 

These conditions are in particular satisfied by the Radau IIA family 
of Runge-Kutta methods [H] . 

The discrete analogue of the variation-of-constants formula 

u(t) = e- tA u + f e- (t ~ T)A g{r) dr 
Jo 

1 Instead of taking r(z) and <?;(z) as rational functions originating from a 
Runge-Kutta method, another suitable choice would be r(z) = e z and qi(z) = 
Jq e*- 1-8 ^ 2 li{8) dd, where ii is the ith Lagrange polynomial corresponding to the 
Gauss nodes Cj. This could be used similarly in the algorithm below. 
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is obtained by solving the recurrence relation ©. With the column 
vector gj = (g(tj + Cih))._ v this becomes 

n-l 

u n = r(-hA) n u + hY, r (-hA) n - l - j q(-hA) gj , n > 1. (9) 
i=o 

3 Discretization of the contour integrals 

We now discretize the Cauchy integral representation 

r(-hA) n q(-hA) = — ( (A + ^l)" 1 r{h\) n q(hX) dX (10) 
2m J r 

along suitable contours r in the resolvent set of —A. The numeri- 
cal integration in (jlOj) is done by applying the trapezoidal rule with 
equidistant steps to a parametrization of a hyperbola [S]. With one 
contour and one set of quadrature points on this contour, we do not 
have a uniformly good approximation for all n = 0, . . . , N, but we 
can instead obtain a uniform approximation locally on a sequence of 
geometrically growing intervals 

I i = [B t - 1 h,B t h), 1>1, (11) 

where the base B > 1 is an integer, e.g., B = 10. For nh G In we 
approximate the contour integrals (J10[) as 

r(-hA) n q(-hA) (12) 

- E wi e \xi e) +A)- l r(hX^rq(hxi e) ), nh G I e , 
k=-K 

with the quadrature points A^, lying on a hyperbola 1^ and with the 

(£) 

corresponding weights w k . The number of quadrature points on i^, 
2K + 1, is chosen independent of I. The contour Fi is chosen as a 
hyperbola given by 

R^T e : 6^j e (e) = M(l-sm(a + iO)) + a (13) 

with an ^-dependent parameter m > 0. The angle a satisfies < 
a < ^ — (f with <p of P)l. and <r is the shift in (|2*|). The weights and 
quadrature points in (|12|) are given by 

= J 7K^fe) , = t^(^) with 9 k = kr , 

where r is a step length parameter that can be chosen independent 
of £ 

The following bound of the necessary number of quadrature points 
is a consequence of the error analysis in Section El 
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Theorem 1 In a quadrature error bounded in norm by e for 

nh £ In is obtained with 

K = Oflog i) . 

T/iis holds for n > c log(l/e), wrai/i some constant c > 0. T/ie required 
number K is independent of I and of n and h < ho with nh < T . For 
o~ < 0, K is also independent of the length T of the time interval. 
K depends on the angle (p, the bound M and the shift a in but 
is otherwise independent of A. 

The approximation is, however, poor for the first few n; cf. also |13j . 

Concerning the choice of parameters we remark that the above 
asymptotic bound for K is obtained with 1/r proportional to log(l/e) 
and with the parameter /j,£ for the contour chosen such that n^B t 'h = 
c\ log(l/e) with c\ independent of I and h, e.g., with c\ = 1/4. Since 
perturbations in the terms of (|12|) can be magnified with r(JiK() n ~ 
e K t nh w [fa K p — ^(l — sin a) + a, the factor c\ should not be cho- 
sen too large. We refer to jHj for an optimized strategy to choose the 
parameters. 



4 The fast algorithm 

We start from the discrete variation-of-constants formula Q for the 
Runge-Kutta approximation with a fixed N. For the expression 
r(hA) N uo we use the discretization of the Cauchy integral like in 
the previous section and in fact similarly to the approach of ^1] for 
computing exp(— tA)uQ. 

The novel algorithm is concerned with the treatment of the inho- 
mogeneity. For a fixed step number N and a given base B we split 
the sum in © into L sums, where L is the smallest integer such that 
N < B L : 

(0) , , (L) 

u N = u K N ' H \-u K N ' 

with uff = hq(—hA)gN-i and 

= h E r(-hA) N ~ l -i q(-hA) 9j 
(N-i-j)helt 

for £ > 1. On inserting the integral representation (|10|) we obtain, 
with m = N - B l for < £ < L - 1 and n L = 0, 

rii—i— 1 

u N= h E 2rij ^ + A)- 1 r{h\) N - 1 -iq{h\)g d\. 
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The integral is discretized with the quadrature formula of Section 03 
we approximate iSS by 17$ given as 

n e _ 1 -l K 

uf =h £ £ wi'W t i)-'K<')»-H 

fe=-Jf 

where 

n^_i-l 

= * £ rrf)"--H g(/( Af) 5 , 

Comparing this formula with Q , we see that yf 1 is the Runge-Kutta 
approximation to the solution at time t = Ui—\h of the linear initial- 
value problem 

y'(t) = \fy{t) + g(t), y(n e h) = 0, (14) 

and hence y£' is computed by Runge-Kutta time-stepping on (fTl|) . 
using © with the scalar hXp in place of the operator —hA. With 
the solutions x k of the linear systems of equations 

(Ai" + A)4«= 9 «, (15) 

we obtain as the linear combination 

<= E W with <£>= „f r^Af)^ 1 . (16) 

There are only (if + 1)L linear systems (|15() to be solved, for = 
0, . . . , K and I < L. (Since the quadrature points lie symmetric with 
respect to the real axis, only the sum of the real parts of half the terms 
in (J16|) needs to be computed when approximating solutions with real 
components.) We recall L — 1 < log B N and K = 0(log(l/e)), where 
e is the accuracy requirement in the discretization of the contour 
integrals. Note that the only approximation made in the computation 
of U N , is the discretization of the contour integrals. 

Because of the poor approximation of the contour integral ()10|) 

for small n, we evaluate uffi + U$ by B direct Runge-Kutta steps 
up to time t = Nh for the initial value problem 

v'(t) +Av{t) =g(t), v((N — B)h) = 0. (17) 
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This requires the solution of another mB linear systems with ma- 
trices of the form (X + A). For small values of B or stringent accu- 
racy requirements, we take B 2 direct Runge-Kutta steps to compute 
vf$ + u$ + ujy . (Asymptotically, we need to take 0(log(l/e)) direct 
steps according to Theorem 1.) 

Finally we sum up the Uj^ to obtain 

U N = U^ + --- + U { N L) (18) 

as the approximation to u^. The fast algorithm thus consists of doing 
the steps (fTl)l - (fTH|) in the given order. 

Remark 1 The algorithm extends to differential equations with a pos- 
itive definite mass matrix M, 

Mu'{t) + Au(t) = g{t), u(0) = n , (19) 

which is transformed to a system u'(t) + Au{t) = g(t) for u(t) = 
M^uit) with A = M^/ 2 AM- 1 / 2 and g(t) = M'^git). Applying 
formally the above algorithm to the transformed system and then 
transforming back yields again (|16|) . where now is the solution of 
the linear system 

{xfM + A) X f=yf, (20) 

((■) 

and y k is the Runge-Kutta approximation at t = ng_\h of the initial 
value problem (|14j) with the untransformed inhomogeneity git). 

Remark 2 We have formulated the algorithm for a constant time step 
size h, but this is not essential. The algorithm is readily extended to 
accommodate variable step sizes, with the same step size sequence 
for all k in (|14[) . chosen adaptively according to the behaviour of the 
inhomogeneity g{t). Adaptivity in space can be used in solving the 
linear systems (|15[1. choosing the spatial mesh according to the be- 

haviour of the right-hand sides y k and the operator A. Note that in 
a hierarchical basis representation, adding a mesh point just corre- 
sponds to adding a scalar differential equation in (|14j) . The details of 
such an adaptive algorithm are beyond the scope of this paper. 
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5 Numerical experiment 

We consider an initial-boundary value problem of the heat equation 
in two space dimensions for u = u(x, t), 

' d t u(x,t) = Au(x,t), x£f2,0<t<T, 

u(x, 0) = 0, x E f2, 

d u u(x,t) = 0, x £ r int , < t < T, 

d v u(x, t) = f3(x, t) - p(u(x, t) - u out ), x £ r out , < t < T, 

on a wire-fence like structure (rectangle of size 10.65 x 12.64 with 
hexagonal holes, each hole with radius 0.8), see Figure ^ Here 
denotes the boundary of the holes, and r out is the boundary of the 
rectangle. In the example we set the heat flux (3 = 5sin 2 (t) on the 
upper and left boundary of the rectangle and j3 = on the lower 
and right boundary, and the convective heat flux to p(u — u ou t), with 
the ambient temperature u ou t = and the coefficient of surface heat 
transfer p = 0.5, cf. the introduction in jTj. Space is discretized using 
linear finite elements on a triangular mesh, with 27346 vertices and 
50368 triangles. Triangulation is done using the tool Triangle |16| . 



Fig. 1. Domain for the heat equation, with isolines of the temperature distribu- 
tion at t = 20. 



The finite element equations are of the form (|19|). where M is the 
standard mass matrix containing the L 2 inner products of the nodal 
basis functions </>j. The stiffness matrix is the sum A = Aq + pMj 
with 

A °\ij = / ^fi^Pjdx > M b\ij = / (fiPjdcr . 
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Fig. 2. Number of solves of linear systems versus step number: direct time- 
stepping (o) and fast algorithm (*). 

The inhomogeneity g(t) is given by 



The algorithm takes into account that g(t) has nonzero entries only 
along the outer boundary r out , so that effectively g(t) is a vector 
whose dimension is the number of degrees of freedom on the outer 
boundary - in this example 776. The differential equations (|14|) need 
to be integrated only for this reduced dimension, since they have no 
coupling between the components. 

We have used the 2- and 3-stage Radau IIA methods (of orders 3 
and 5, respectively) for time discretization in our numerical experi- 
ments. 

In the fast algorithm we set B = 5 and K = 15 and, from the 
experience of [9|ll3| . we choose the angle in the hyperbola as a = vr/4, 
the parameter fi£ = 3/(hB e ) and the parameter r = 5/K. This choice 
of parameters leads to a deviation of the order 10 -6 from the Runge- 
Kutta approximation at time t = 20. 

The two-dimensional example above is still small enough that a 
direct solution of the linear systems using sparse solvers is reasonable. 
A direct implementation of the m-stage Radau IIA method (cf. [SJ) 
requires only m sparse LU factorizations, computed at the beginning 
of the integration, followed by mN substitutions. On the other hand, 
for the algorithm presented here we need to solve (K + 1)(L — 1) 
linear systems with matrices AM + A for as many different values 
of A, and the mB linear systems for the B direct steps. Especially 
with a diagonal, lumped mass matrix M = DD T , this work can be 
reduced by a similarity transform taking D~ 1 AD~ T to tridiagonal 
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(or Hessenberg) form T, but exploiting sparsity here becomes an is- 
sue; see |Hll2|. The resulting linear systems with XI + T are then 
inexpensive to solve. Even without using such a transform, the fast 
algorithm eventually overtakes the standard algorithm for sufficiently 
large step numbers N, in the present example for N ~ 1000. Much 
earlier and larger relative gains arise when iterative solvers are used 
for the linear systems in both algorithms, as is clear from the linear 
systems count in Figure 12 



6 Error analysis 

Our analysis relies on the good behaviour of the trapezoidal rule 
for certain holomorphic integrands [Hl ll7lll8j . Following the ideas in 
[5j, we consider the continuation of the parametrization (|13|) to the 
conformal mapping 

7(tt>) = fj, (1 — sin(o + iw)). (21) 

(For ease of presentation we set a = in (J2J.) This conformal map- 
ping transforms each horizontal straight line 

Imto = y, — d < y < d, 

with < a — d < a + d < ^, into the left branch of the hyperbola 

A C- ( ReA_/i V ( ImA V - i 
\lism(a-y)J \ficos(a-y)J 

i.e., the left branch of the hyperbola with center at (//, 0), foci at 
(0,0), (2/x, 0) and with asymptotes forming angles ±[vr/2 — (a — y)] 
with the real axis. Therefore, 7 transforms the horizontal strip 

D d = {w£C: \lmw\ < d} 

into the region I? = j(Dd) limited by the left branches corresponding 
to y = ±d. To indicate the dependence on the parameter fi of ()21|) . 
we write Q = fl^. We note that A G i?^ if and only if hX £ for 
any h > 0, so that 

hQ^ = {2^. 

Because of (J2J), henceforth we will assume that a > and d > 
satisfy < a — d < a + d < \ — ip. Under these conditions, all the 
hyperbolas we are considering lie outside the spectrum of —A. 
After parametrizing (|10j) via 7, we get 

r+00 

r(-hA) n q(-hA) = / G h Jx)dx, 
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where Gh, n {w) is given, for w 6 Dd, by 

G h>n (w) = ^- (7H + A) r(fry(«0) n T 'H- (22) 

For an integrable mapping G : R — ► X, iT > 1 and r > 0, set 

/+00 
G(x)dx-r G(fcr), (23) 
■°° fe=-K 

i.e., E t k(G) stands for the quadrature error of the truncated trape- 
zoidal rule for the integral of G. Our goal is precisely to estimate 
E Tj K{Gh,n)- To this end we first consider the behaviour of G^n on 
Dd- We need the following lemma whose elementary proof is omitted. 

Lemma 1 Let r{z) be a rational function with r(0) = 1, r'(0) = 1 
which satisfies the L-stability condition f%|). Then, there exist p > 
and b > such that 

e 2S 

\r(z)\ < 1 + f° r z(^Q 5 with < 5 < p. (24) 

Now, from the sectorial condition (J2J) on A and Lemma ^ with 
5 = hp < p, we obtain 

2fihn 

\\G h , n (x + iy)\\<C — — (25) 

(1 + onp(coshx — sin (a — y))) 

for x E R and \y\ < d (recall that < a — d < a + d < ^ — </?), where 
Cq is the constant given by 



_ M l + sm(a + d) „ , , 

C = 7r\hi ~t — n\ max||g(z) 



Finally, the above bound (|25|). the elementary inequality 

l + c-s> (l-s)(l + c), c,s>0, 
and the bound 1 for the sine yield, for |y| < d and i = n/i, 



-n 



(J ^2fj,t . but 
\\GU* + W)\\< (1 .° M/B) . (l + ^coshxj . (26) 

Next, to estimate E T: K{Gh^ n ), we are going to use an approach 
similar to the one in [Hlll7lll8j . We denote by S(Dd,X) the class 
formed by all the continuous mappings G : Dd — > X (for a complex 
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Banach space X, here a space of matrices) holomorphic on the interior 
of the strip D d , and satisfying the following two conditions: 



j \\G{x + iy)\\dy -> 0, as \x\ -> +00, (27) 

J-d 



N(G,D d ):= {\\G(x + id)\\ + \\G(x - id)\\} dx < +00. (28) 

J —00 

Given G G S(D ( i,X), it turns out, assuming that G has a fast 
decay at 00, that E Tt K{G) becomes very small as K —* +00 if r is 
properly tuned (see [5l ll7lfT5] for various situations). In Theorem [2] 
we assume that G exhibits the kind of decay of G/j jn in (|26|) and 
this theorem will directly provide the estimate for E T ^K{Gh,n) we are 
looking for. 

Theorem 2 Assume that G £ S(D d , X) for some d > 0, and that 
there exist C, a > and n > 1 suc/i £/ia£ 

||G(x)|| < C(l+ -coshxV™, (29) 



T/ien, /or r > 0, K > 1, there holds 
N(G,D d ) 



\E t ,k(G)\\ < 



e 2-Kd/r _ \ 



+C (^(a)e- acosHKT)/2 + ( 



a \ -(n-l) 

1 H cosh ivr 

n 



with 4>(a) =2 + | log(l — e 



-a/2} 



Notice that <fi is decreasing, (p(a) — ► 2 as a — ► +00 and 0(a) ~ 
I log a| as a — > 0+. 

Proof Set 

/+00 00 
G(x)dx-r G(frr), r > 0. 

-°° fc=-oo 

For fixed A' > 1, it is clear that 

\\E TjK (G)\\<\\E Tj00 (G)\\+t W G ( kT )l 

\k\>K+l 

On the one hand, by Theorem 4.1 in (see also [E])> we have 

N(G, D d ) 



l^r,oo(G)|| < 



e 2ird/r _ ^ ' 
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On the other hand, 



r \\G(kr)\\ < 2Ct ^ (l + -coshfcr 

|fe|>if+l k=K+l " 



<2C 1 + -coshx dx. 

Jkt v n J 

The proof of the theorem is now completed by applying the following 
lemma. □ 

Lemma 2 For R > 0, a > and n > 1 there holds 

(l + ^coeh < ^(a)e~ acoshiJ/2 +(l + ^coshi?) _(ri_1) . 

Proof The change of variables u = cosh x shows that 

a \-n f + °° ( a \~ n du 

1 H — cosh x) dx = I 1 H — u , 



R 



n ' JcoshR ^ n J vV - 1 ' 



Set (5 = max{ cosh R, n/a}. Then, from the estimates in jS] and the 
elementary inequality 

(1 + y/n)~ n < e~ y/2 , for < y < n (30) 

it turns out that 

l + -u < / e au/ - 



coshi? x n ' Vu 2 - 1 JcoshR Vu 2 - 1 

/•+oo 

< / e - acoshx/2 dx 

< (/.(a)e- acoshfi/2 . 

Moreover, 

h ^ ' a \ ~ n du 
1 + -w 



/J 



n J ^Ju 2 - 1 
^(l + ^coshi?)^/ 



(n-i) /■+« / a \ -i 
l + -u] 



n ) yju 2 - 1 



Now, since /? > max{l, n/a}, the result follows from the observation 
that for both n/a > 1 and n/a < 1 we have 

/•+ 00 / a \-i /• +00 , 1 . i dv 

/ l + -u < / (l + vy 1 — r = 1- □ 

imax{l,n/a} V n ' V U 2 - 1 A Vf 2 -1 
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We apply Theorem |21 to Gh, n - First of all, notice that by (|2fi[) it is 
clear that Gh, n satisfies l|27[). Moreover, by Lemma El we have 

4Cn e 2 ^ 

N(G hn , D d ) < , , — . (31) 



^t)e-^ + (l + ^)~ (n - 1) ). 



and conclude that £ S(Dd,X). Then, in view of (|26j) and (|31jl. 
Theorem |2] yields directly 

ii ^ mi / 4C o e 2Mi / <KM) e"^ /2 + (1 + biit/n)-^- 1 ) 
\\ErMG h ,n)\\ < (1 _ M/n)n ^ ^57731 

A simplified version of this estimate is obtained by using the elemen- 
tary inequalities (|3(Jj) and 

(1 - y/n)' n < e 2y , for < y < n/2, 

(p(y) < 3, for y>l. 

Setting 

C = 20 C , a = 2 + §6, ai = 2 + 26, a 2 = 

with 6 of Lemma ^ as before, we can summarize the final result in 
the following theorem. 

Theorem 3 The quadrature error (JUJ) for G h>n of {H§ wii/i flUJ 
satisfies, for t = nh and if n/2 > but > 1, 

( p ao)J,t 
e 2«d/T_ 1 +e 

+ e^(l + ^cosh(Kr))^ 1) y 

The first term in the error bound becomes 0(e) if r is chosen 
so small that aofit — 2nd/ t < loge, which requires an asymptotic 
proportionality 

1 , 1 

- ~ log - + fit. 



16 



M. Lopez-Fernandez, C. Lubich, C. Palencia, and A. Schadle 



For fi chosen such that 

ci, 1 ,1 

— log - < < a log - 

Be e 

with an arbitrary positive constant c\ and with B > 1, we obtain 
that the second term is 0(e) if a\ — ci2 cosh(i^r) < —B/c\, i.e., with 

cosh(Kr) = C2 

for a sufficiently large constant ci. With the above choice of r, this 
yields 

K ~ log - . 

£ 

The third term then becomes smaller than e for 

n> c log — 

with a sufficiently large constant c. Taken together, these estimates 
prove Theorem^ 
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